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EIGENVECTOR DYNAMICS: GENERAL THEORY AND SOME 

APPLICATIONS 

ROMAIN ALLEZ^'^ AND JEAN-PHILIPPE BOUCHAUD^ 



Abstract. We propose a general framework to study the stability of the subspace 

spanned by P consecutive eigenvectors of a generic symmetric matrix Ho , when a small 

perturbation is added. This problem is relevant in various contexts, including quantum 

dissipation (Ho is then the Hamiltonian) and financial risk control (in which case Ho is 

the assets return covariance matrix). We argue that the problem can be formulated in 

' ' terms of the singular values of an overlap matrix, that allows one to define an overlap 

^ distance. We specialize our results for the case of a Gaussian Orthogonal Hq, for which 

(^ the full spectrum of singular values can be explicitly computed. We also consider the case 

O when Ho is a covariance matrix and illustrate the usefulness of our results using financial 

data. The special case where the top eigenvalue is much larger than all the other ones can 

be investigated in full detail. In particular, the dynamics of the angle made by the top 

eigenvector and its true direction defines an interesting new class of random processes. 
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1. Introduction 

Random Matrix Theory (RMT) is extraordinary powerful at describing the eigenvalues 
statistics of large random, or pseudo-random, matrices |27t [2l [3ll6]. Eigenvalue densities, 
two-point correlation functions, level spacing distributions, etc. can be characterized with 
exquisite details. The "dynamics" of these eigenvalues, i.e. the way these eigenvalues evolve 
when the initial matrix Hq is perturbed by some small matrix eP, is also well understood 
[7] . The knowledge of the corresponding eigenvectors is comparatively much poorer (but 
see [22] )• One reason is that many RMT results concern rotationally invariant matrix 
ensembles, such that by definition the statistics of eigenvectors is featureless. Still, as 
we will show below, some interesting results can be derived for the dynamics of these 
eigenvectors. Let us give two examples for which this question is highly relevant. 

One problem where the evolution of eigenvectors is important is Quantum Dissipation 
[50] (see also the related recent strand of the literature on Quantum "Fidelity" [8]). As 
the parameters of the Hamiltonian Hj = Hq + sPt of a system evolve with time t, the 
average energy changes as well. One term corresponds to the average (reversible) change 
of the Hamiltonian that leads to a shift of the energy levels (the eigenvalues). But if the 
external perturbation is not infinitely slow, some transitions between energy levels will 
take place, leading to a dissipative (irreversible) term in the evolution equation of the 
average energy of the system. The adiabaticity condition which ensures that no transition 
takes place amounts to comparing the speed of change of the perturbation ePt with a 
quantity proportional to the typical spacing between energy levels. For systems involving 
a very large number N of degrees of freedom, the average level spacing of the N x N 
Hamiltonian H goes to zero as N~^. For N — )• oo, any finite speed of change therefore 
corresponds to the "fast" limit, where a large number of transitions between states is 
expected. In fact, if the quantum system is in state \(j)^) at time t = 0, that corresponds 
to the ith eigenvector of Hq, the probability to jump to the jth eigenvector of Hi, \(l)]), 
at time t = 1 is given by |(0i|</'i')P, where we use the bra-ket notation for vectors and 
scalar products. The way energy is absorbed by the system will therefore be determined by 
the perturbation-induced distortion of the eigenvectors. More precisely, if \(j)'^) is different 
from \(l)j), some transitions must take place in the non-adiabatic limit, that involve all the 
states j which have a significant overlap with the initial state. 

Another very relevant situation is Quantitative Finance, where the covariance matrix 
C between the returns of A^ assets (for example stocks) plays a major role in risk control 
and portfolio construction [S]. More precisely, the risk of a portfolio that invests Wa in 
asset a is given by 7^^ = X^Q.ait'QCa/^tt'/j. Constructing low risk portfolios requires the 
knowledge of the n largest eigenvalues of C (n is often chosen empirically, keeping only 
the statistically meaningful eigenvalues that lie outside the Marchenko-Pastur sea, see [13] 
for details), Ai ^ ... ^ A„ and their corresponding eigenvectors \(pi), ■ ■ ■ , \4>n)- The top 
eigenvalues and eigenvectors represent the most risky directions in a financial context. A 
portfolio such that the vector of weights \w) has zero overlap with the first n eigenvectors 
of C has a risk that is bounded from above by A,i+i. The problem with this idea is that 
it relies on the assumption that the covariance matrix C is perfectly known and constant 
in time. The observation of a sufficiently long time series of past returns would thus allow 
one, in such a stable world, to determine C and to immunize the portfolio against risky 
investment modes. 

Unfortunately, this idea is thwarted by two (inter-related) predicaments: a) time series 
are always of finite length, and lead to substantial "noise" in empirical estimates of C [9] 
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and b) the world is clearly not stationary and there is no guarantee that the covariance ma- 
trix corresponding to the pre-crisis period 2000-2007 is the same as the one corresponding 
to the period 2008-2011. For one thing, some companies disappear and others are created 
in the course of time. But even restricted to companies that exist throughout the whole 
period, it is by no means granted that the correlation between stock returns do not evolve 
in time. This is why it is common practice in the financial industry to restrict the period 
used to determine the covariance matrix to windows of a few years into the recent past. 
This leads to the measurement noise problem alluded above. Now, if the "future" large 
eigenvectors do not coincide with the past ones, a supposedly low risk portfolio will in 
fact be exposed to large risks directions in the future. Denoting as \(f)^) the past eigenvec- 
tors and \(j)^) the future ones, the total risk of the portfolio \w) = |(/)^) can be defined as 
X]?=i ^]{4']\4'^)'^- Therefore, as for the quantum dissipation problem, the statistics of the 
overlaps {(j^^Acj)^) is a crucial piece of information. 

In practice, one computes the empirical covariance matrix E using past stock returns, 
which is defined as: 

1 ^ 
t=i 

where T is the length of the period on which the measurement is done and r\ is the return 
of stock number i on day t. If the true covariance matrix C exists and is stable in this 
period, the empirical matrix E can be seen as a perturbation of C, since one can write 
E = C + iS where £ \s a. matrix whose elements are of order 1/vT (by the central limit 
theorem) to be considered small as T is usually quite large. In this sense, the problem falls 
in the more general context introduced above. 

The paper is organized as follows. In the next section[2| we introduce the main statistical 
tools and objects studied in diff^erent contexts in the following sections and we also briefly 
recall standard perturbation theory. In the next two sections, we turn to two explicit 
illustrations, first in the context of matrices in the Gaussian Orthogonal ensemble (GOE), 
and then in the context of covariance matrices. More precisely, in section [3j we study 
the stability of the eigenvectors for a matrix Hq in the GOE by computing the "overlap 
distance" between the perturbed space and the non-perturbed space in the limit of large 
matrices Hq when the perturbation matrix P is also in the GOE. Furthermore, we are able 
to compute the full spectrum of the overlap matrix in this limit, which gives a precise idea 
of the perturbation induced distortion for the eigenvectors. In section |4j we go through 
the same steps in the context of covariance matrices. We study the link between the 
population eigenvectors (the eigenvectors of the true covariance matrix) and the sample 
eigenvectors (the eigenvectors of the empirical covariance matrix) . Then, in section ^ we 
analyse more precisely the case of a population covariance matrix with an isolated top 
eigenvalue much larger than the other ones. We measure the empirical covariance matrix 
with an exponential moving average estimator and characterize the temporal evolution of 
the angle made by the top eigenvector and its true direction which defines an interesting 
new class of random processes. Finally, in section |6j we apply our ideas to the analysis of 
financial market correlations. Our purpose here is to study whether correlations between 
stock returns evolve or not. In particular, is there a constant in time correlation matrix 
(population correlation matrix)? Do the economical sectors (eigenvectors of the correlation 
matrix) evolve or not ? We find that there is indeed a genuine evolution of the correlation 
matrix of stocks returns for different markets in the U.S, in Europe and in Japan, a result 
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that confirms recent studies (see e.g. [lOl HU 12] )• We also give a partial description of 
this temporal evolution. 

2. Perturbation theory and Statistical tools 

In the whole paper, we will mainly be interested in the eigenvectors of a matrix Hi that 
can be written as 

Hi = Ho + eP 

where Hq and P are two N x N symmetric matrices and e a small (positive) parameter. 
The matrix Hq is the true signal which is perturbed by the adding of the small term eP. 
The matrix Hi will be referred as the perturbed matrix. The eigenvalues of the matrix 
Hj, i = 0, 1 will be denoted as X\ '^ X\ '^ ... ^ A'^y and the corresponding eigenvectors 

l'Al>,...,l05v)- 

Our aim is to describe the relation between the perturbed eigenvectors \<j)j) and the 
non-perturbed eigenvectors {(p^) when the parameter e tends to 0. 

When trying to follow the evolution of a given eigenvector \(j)i) when the small pertur- 
bation eP is added, one immediately faces a problem when eigenvalues "collide" . It is well 
known that true collisions (degeneracies) are non generic; the collisions are in fact avoided 
and levels do not cross. However, upon the pseudo-collision of Aj and Aj+i (say), the eigen- 
vectors \(j)i) and |i;^i-i-i) strongly hybridize (this phenomenon was observed for example in 
|28t Fig. 1]). The eigenvector \(p^) will in fact hybridize with all the eigenvectors associ- 



ated to eigenvalues Xpj 7^ i that are close to A^ and will therefore be strongly affected 
by the adding of the perturbation if the level spacing is too small. Indeed, using standard 
perturbation theory to second order in e, the perturbed eigenvectors can be expressed in 
terms of the initial eigenvectors, for small e, as: 

(2,1) 




+-'^E 



V^ ^ji^ii i'^iit^ij \ I ,0\ 



where Pij = ((/>^|P|(/)^). The denominators A? — A^ remind us that eigenvectors are strongly 
affected by eigenvalue pseudo-collisions, as alluded to above. We mention in passing that 
perturbation theory to second order in e for the eigenvalues gives 

X} = X^ + eP,, + e^Y.^Jk_^ (2.2) 



It is important to mention at this point that equations (2.1) and (2.2) are a priori only valid 
in the regime where the entries of the perturbation matrix eP are small compared to the 
level spacing of the non-perturbed matrix Hq. This condition ensures that the asymptotic 



correction terms appearing in (2.1) and (2.2) are small compared to the leading term of 



order 1 corresponding to the non-perturbed system. 

The idea is then to study the stability of a whole subspace Vq spanned by 2p -|- 1 several 
consecutive eigenvalues: {\(p^ ), ■ ■ ■ \(t>\), • • ■ , I'/'fc+p)}- Motivated by the above examples, 
we ask the following question: how should one choose q ^ p such that the subspace Vi of 
dimension 2q + l spanned by the set {{(pk-q)^ • • • 1*^1)' • • • ' \^k+q)} ^^^ ^ significant overlap 
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with the initial subspace? In order to answer this question, we consider the (2q+l) x (2p+l) 
rectangular matrix of overlaps G with entries: 

G,r-= {^}\^°j) ■ (2.3) 

The {2p + 1) non zero singular values 1 ^ si ^ S2 ^ ... ^ S2p+i ^ of G give full 
information about the overlap between the two spaces. For example, the largest singular 
value si indicates that there is a certain linear combination of the {2q + 1) perturbed 
eigenvectors that has a scalar product si with a certain linear combination of the {2p + 1) 
unperturbed eigenvectors. If S2p+i = 1, then the initial subspace is entirely spanned by 
the perturbed subspace. If on the contrary si <C 1, it means that the initial and perturbed 
eigenspace are nearly orthogonal to one another since even the largest possible overlap 
between any linear combination of the original and perturbed eigenvectors is very small. 
A good measure of what can be called an overlap distance D(Vo,Vi) between the two 
spaces Vq and Vi is provided by the average of the logarithm of the singular values: 

^(^o,V^i):=-§^, (2.4) 

but alternative measures, such as 1 — Yli^i/C^P + l)i ^^^ be considered as well. Since 
the singular values s are obtained as the square-root of the eigenvalues of the matrix 
G^G, one has D = — Indet G^G/2P, where we henceforth introduce for convenience 
the notations P = 2p + 1, Q = 2q + 1. The overlap distance D was originally studied 
for P = (5 in |lj, see e.g. [151 128], where a fundamental effect observed in many body 
systems, called the Anderson Orthogonality catastrophe (AOC) is introduced. Anderson 
[1] addressed the ground state of a finite system consisting of P noninteracting electrons. 
Upon the introduction of a finite rank perturbation matrix eP, this ground state gets 
modified. It is then shown that the overlap between the original and the modified P- 
electron ground state, which is in fact exactly given by our overlap distance D{Vo,Vi) 
between the two subspaces Vq and Vi (with P = Q), is proportional to a negative power 
of P, and vanishes in the thermodynamic P — )• +00 limit, hence the catastrophe. We will 
see that our idea of introducing a rectangular Q x P overlap matrix G enables to avoid 
this orthogonality catastrophe. Our objects introduced here will also allow us to revisit the 
AOC in the case of square matrices G showing that it occurs for the random matrix model 
studied in section |3] (AOC for this RM model is also studied in [28]). In [13 [28], the AOC 
is also investigated through random matrix models as in our paper. The main difference 
with [15j is that we consider here full rank perturbation instead of a localized perturbation 
of rank 1, for which one can do explicit computations (and so treat the non-perturbative 
regime) . 

As an interesting benchmark, consider the case when two subspaces Wq and Wi respec- 
tively of dimensions P and Q are constructed using randomly chosen orthonormal vectors 
in a space of dimension A^. In this case, one expects accidental overlaps, such that the Si 
are in fact non zero, and therefore D{Wo,Wi) is finite. This distance can be calculated 
exactly using Random Matrix Theory tools in the limit N,P,Q ^- 00, with a = P/N and 
(3 = Q/N held fixed. The result is [II]: 



Jo pTTS[i — S ) 

where 7-1- = a -|- /3 — 2a/3 it 2y^al3{l — a)(l — /3). In other words, in that limit, the full 
density of singular values is known; all singular values are within the interval [a/t^, v/7+]. 
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This provides a benchmark to test whether the two eigenspaces are accidentally close 
{D « Drmt), or if they are genuinely similar (D <C Drmt)- 

Endowed with the above formalism, we can now proceed to compute -D(Vo, Vi) in the 



case where the perturbation is small. Indeed equation (2.1 ) allows us to obtain the overlap 
matrix G. Keeping only the relevant terms to order e^, we findjj 

2 



Gij 



1 ^2 S^^i ( \0-\0 ] if « - j 



2 ^l^i VaO-AO , , 

"" - ^2.5) 



Using (2.5), we can also compute the matrix G^G to second order in e, we obtain for 
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i / j: 



and, for i = j: 



{G^G)u = l-e^ Y^ U^ . (2.7) 

It is then easy to derive the central result of our study: to second order in e, the distance 
L>(Vb) ^i) between the initial and perturbed eigenspaces is: 

i=k—pj^{k—q,...,k+q} \ 1 * / 

The matrices G and G^G are both close to the identity matrix as they should. Let us 
define the matrix S by 

S = ^ fl - GtG 

One can note that the matrix S is positive definite and that its matrix elements are of 
order 1 when e goes to 0. 

3. Eigenvector stability in the GOE ensemble 

We will now define a random matrix model for which we will apply the results of the 
previous section. Let Hq be a random matrix of the Gaussian Orthogonal Ensemble (GOE), 
i.e. a matrix of size N x N with gaussian entries randomly chosen with the probability 
measure on the space of real symmetric matrices 

N 

This definition implies that the matrix Ho is symmetric with independent Gaussian entries 
above the diagonal with variance a'^/N on the diagonal and a'^/2N off diagonal. 

The perturbation matrix is similarly defined as a random matrix of the GOE, indepen- 
dent of Hq with the same variance profile for the entries. 

We then define the perturbed matrix Hi as before: 

Hi = Ho + eP. (3.1) 



P{dHo) = eM-:^HHi)) dHo . 



see [29| for similar calculations. 
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It is very well known that the density of Ho-eigenvalues PnW '■= l/-^l^j=i ^A^ tends in 



the large A^ limit to the Wigner semi-circle law 

p{dX) = -'-\/4(T2-A2dA. (3.2) 

ZTT 

For simplicity, we take cr^ = 2 in the following. 

Remark 3.1. Here the choice of a GOE random matrix for Ho is made to get an explicit 
expression for the density of states in the limit of large matrices. But our theory developed 
in the following would apply for sequences of matrices (Ho(A^))Ar such that the density of 
states converges to a general (not necessarily the semi-circle density) continuous density 
p(A)dA. Moreover, the sequence (Ho(A^))Af can he supposed deterministic or random. The 
matrix Hq can he seen as the true signal to which a small noisy perturbation Pj is added. 

In this whole current section, ^~^ denotes an averaging over the random matrix P q 

3.1. Distance between subspaces of perturbed and non-perturbed eigenvectors. 

We consider the subspace Vq of initial eigenvectors corresponding to all the eigenvalues A 
contained in a certain finite interval [a, b] included in the Wigner sea [—2, 2]. We want to 
compute the mean overlap distance D{Vo,Vi) between Vq and the subspace Vi spanned by 
the perturbed eigenvectors of Hi, corresponding to all eigenvalues contained in [a — 6, b+6], 
where 5 is a positive parameter. 



Using formula (2.8), which is valid if the entries of the perturbation matrix eF (of 
order eN~^i'^) are much smaller than the mean level spacing of the matrix Hq, of order 

{N p{\))~^., we can write for e <^ N^^l"^: 

D{Vo,V,) = ^ ^ J: -,^. (3.3) 



It is easily seen that Eq. (3.3) becomes, in the large N limit: 



D(V.V,) = -^ fax f dVftM^. (3.4) 

2 /^ p{X)dX J a J[-2;2]\[a-5;b+5] (^ - A')^ 



where p is the Wigner semicircle density (3.2) 



Formula (3.4) is a priori only rigorously valid in the perturbative regime where e <C 
_/V^i/2, \Ye argue that in fact it remains valid in a wider regime where e ^ 1. Indeed 
although perturbation theory for the eigenvectors fails for Hq eigenvalues that are at 
distance of order of the mean level spacing of Hq, it remains valid in the limit e <C 1 for 
eigenvalues at distance large compared to the order of the perturbation entries £N~^''^ 
and in particular for two eigenvalues lying respectively in the two well separated intervals 
[a; b] and [a — (^; 6 + 5] for which this distance is larger than 6 (which indeed is ^ eN~^''^). 



We see that every terms appearing in (2.1) corresponding to overlap between eigenvectors 



associated to eigenvalues that are at distance smaller than 6 disappear in formulas (2.8) 



(and also in (2.6), (2.7)). Therefore, we expect (3.4), as well as our results below, to remain 
valid in the regime N~'^''^ <C e <C 1, provided the computed distance D{Vo,Vi) itself 
remains much smalleilj than unity. We checked formula (3.4) using numerical simulations. 



with very good agreement for different values of a, b, 6, N, e. In those numerical tests we 



There is no need in averaging over the random matrix Ho for the following results to be valid. 

For this condition to be valid, 5 has to be fixed independent of e, or at least such that e^| ln((5)| <C 1. 
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chose the parameters A^, e, 5 so as to approach the regime N~^''^ <^ £ <^ 1 (for example, 
iV = 4000, £ = 0.1, 5 = 0.5). _ 

We will now write D{a, b; d, e) instead of D{Vo, Vi). 

It is interesting to study the above expression in the double limit (5—7-0 and A = 6— a — )■ 
0. One finds: 

l-n, » . > ( "'""a^^" if««A«l, 

-,D(a,a + A:S..)^^^- if A « * « 1. '^'■'^ 

This last expression shows that when the width A of interval [a, b] tends to zero, the 
corresponding eigenvectors are scattered in a region of width 5 much larger than A itself 
as soon as e ^ vA. It also shows that for fixed A, the distance D diverges logarithmically 
when (5 — )■ 0. This is a consequence of the pseudo-collisions that occur between eigenvalues 
close to the boundaries of the interval [a,b]. When 5 > 0, these pseudo-collisions are 
avoided and D remains finite. 



When 6 = 0, we can do a more precise analysis of the right hand side of (3.3). One can 



show, for large N, that the following result holds at least in the regime e <^ N ^''^: 

\D{a, b;5 = 0,e)^\nN ^^"'\ ^ ^^^^ + A{a, b) (3.6) 

where A(a, b) is a constant independent of N that can be explicitly computed, and involves 
the well known two-point function g{r) that describes the level-level correlations in the 
GOE (see Appendix A for the details of this computation). The IniV term can be guessed 
from the logarithmic behavior of D when 5 — )• 0, since one indeed expects the divergence 



to be cut-off when 6 becomes of the order of the level spacing, i.e. 5 ~ {Np)~ . Eq. (3.6) 
is precisely the Anderson orthogonality catastrophe as first introduced in p!] in the case 
of finite rank perturbation matrices. We recover here exactly the result of [25j (see their 
Eq. (31)) by taking a = -2, 5 = in our Eq. ^^. 



As a side remark, we note that Eq. (3.5) predicts that a fraction oc 5 ^ of the original 



eigenspace gets shoved away at distances larger than 5 (in eigenvalue space) . In the context 
of the non adiabatic evolution of a quantum system [30] , this implies that the energy of the 
system makes jump with a power-law distribution of sizes, such that all moments of order 
q ^ 1 diverge. This means that under an extreme non-adiabatic process, the energy is not 
diffusive but rather performs a "Cauchy flight" (i.e. a Levy flight with a tail exponent 
equal to 2), see p9.J. When the perturbation varies on a flnite time r, one expects the tails 
of the jump process to be truncated beyond 5 ~ h/r [30] . 



3.2. Full distribution of the singular values of the overlap matrix. To order e^, the 
distance D computed in the previous subsection is proportional to the mean position of the 
singular values. One can actually be much more precise and compute, for P, Q — )• cxd, the 
full distribution of all singular values, giving an indication of their scatter around the mean 
position (s). The computation of the density of states (DOS) can be straightforwardly 
performed using free random matrices techniques. 



The authors of [28] expect deviations in (3.61 when the parameter x := £\/iV (which has to be ^ 1 



for ( |3.6| l to be fully valid) is increased. However their numerical results (presented in Fig. 2 of [2H|) show 
that the discrepancies are only noticeable for x close to 1. In addition, the authors of [28] explain that the 



failure of (3.61 for not small enough x is due to the first-order perturbation theory estimate that breaks 
down when used for levels in the vicinity of the edges a, b of the initial interval. This problem was avoided 
previously by the use of rectangular matrices with Q > P and the introduction of the 5 margin at the 
edges a and b. 
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Using Eq. (2.6) and (2.7), it is easy to compute the entries of the matrix 5] in the 
perturbative regime £\/N <^ 1: [j 

^-= ^ (X xfx X) ^'-'^ 

oaii, i,^ X '■-^^ ~ Xi,)[Xj - Xi) 

l^{k-q\...;k+q} 

Denote, for each I ^ {k — q] . . .■,k + q] by V£ the random Gaussian vectors of M^ 

_ / Pe,k-p Pe,k-p+i Pe,k+p 

\Xk-p — Xi' Xk-p+i — A^' ' Xk+p — Xe^ 

It is easily seen that in fact (changing to the equivalent notation for the summation on 
i in term of a, b) 

t.Xii^[a-5;b+5] 

This matrix v^vj is clearly the matrix of a projector on v^ and has only one non-zero 
eigenvalue which is equal to 

.(A,)H|v.lli^ E (jf^S 

je{k-p;...;k+p} ^ ^ "^ 

which can be approximated in the limit of large matrices (P — >■ oo) by 

a(A,)^ f'dX- ^^^^ 



(A-A,)2- 
The resolvent Z^{z) = ;ptr((z — v^vj)^^) of the matrix v^vj is equal to: n 

The Blue function, which by definition is the functional inverse of the resolvent B {Z (z)) 
z, can be computed to first order in 1/P: 



z Pl-a{Xi)z 
Finally, the Red function, defined as R (z) = B (z) — -, is given by: 

''^'> Pl-a{X,)z 

The trick, coming from the theory of free matrices, is to use the additive property 
of the Red function (also called R-transform) for the asymptotically free matrices v^vj. 
Essentially, the R-transform of the matrix 5] can be computed as the sum of the R- 
transforms of the matrices v^ vj : 



e<^{k-q;...;k+q} e^{k-q;...;k+q} 



a{Xi 



-q;-;k+q} e^{k-q;...;k+q} ^ "^ 



We skip the subscript on the eigenvalues A^s. 

Resolvents are usually denoted by the letter G, but we do not want to confuse the reader with the 
overlap matrix G of which we compute the singular value spectrum. 



10 



ROMAIN ALLEZi'2 AND JEAN-PHILIPPE BOUCHAUDi 



Finally, the Blue function of X) is: 

Biz) 



1 1 

z ^ P ^^ 1 

e'^{k-q;...;k+q} 



E 



a{Xi 



which can be approximated in the limit of large P as: 



B{z) 



1 1 



dX 



a{Xi)z 



p{X)a{X) 



-2;2]\[a-<5;6+5] 1 " a{X)z 



(3.8) 



where we note here and below N^ := f p{X)dX. Rewriting equation (3.8) in terms of the 
resolvent gives our central result: 



1 1 



dX 



-2;2]\[a-S;b+5] 



P(A)^(A) 
1 - ct{X)Z{z) ■ 



(3.9) 



Equation (3.9) characterizes the density of states of the matrix Xl in the limit of large 



dimension. We ran numerical simulations to test the validity of Eq. (3.9) in the regime 
1/yN <C e <C 1, see Fig. Ill The agreement is excellent. It would be interesting to run this 
numerical test for very large values of N (here we took N = 4000) so as to fully reach the 
regime 1/yN ^ e ^ 1. However, this becomes numerically demanding, and we leave this 
study for future work. 

We now want to extract the qualitative informations about the distribution of all sin- 
gular values of the matrix G from this equation. In particular, we will show in the next 
subsection that the density of singular values has a compact support for which we char- 
acterize the left and right edges. We also study the shape of this distribution in the two 
asymptotic regimes A ^ 5 and 5 <C A <C 1. 

3.3. Qualitative properties of the spectrum of S. 

Right and Left edges. The relation between the resolvent Z and the density of states r{s) 
of the matrix 5] is limt^;_s.o ^^(s — ioj) = 7rr(s), s G M. Note that one should not confuse 
the density of states p{X) of the original matrix Hq with the density of eigenvalues r{s) of 



From equation (3.9), we can derive a system of equations for the real {g{s)) and imagi- 

p{x)a{x){l - a{x)g{s)) 



nary (r(s)) parts of Z{s), for a; — )• 0: 



9{s) 



+ 



1 



dx- 



(7(s)2+^2^(s)2 ATb J (l-a(x)(7(s))2 + a(x)27r2r(s)2 

[~2;2]\[a-S;h+5] 



(3.10) 



= r{s) 



-1 



g[s)'^ + 7r^r(s) 



+ 



1 



dx- 



p{x)a{xY 



\ 



-2;2]\[a-(5;fe+5] 



1 — (T{x)g{s)Y + cr(x)27r2r(s)2 



(3.11) 



The second equation (3.11) always admits the solution r = 0. Plugging r = into the 
first equation gives: 



1 1 

W)^Nb 



dx 



[-2;2]\[a-5;b+&] 



1 -a{x)g{s) 



(3.12) 
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Figure 1. The histogram represents a numerical simulation of the density of 
states of the matrix S (computed with 15 independent samples). The red curve 
is the theoretical corresponding density for r{s) obtained by solving numerically 
( 3.9 ). For this figure, we chose a = 0,b — 0.5, 6 — 0.5. We chose the parameters N 



and e so as to approach the "less perturbative" regime where 1/yN ^ e ^ 1 for 
this figure as iV = 4000 and e = 0.1. 



Equation (3.12 ) implies the asymptotic relation g{s) 
values of s correspond to small values of g[s). Set 



mo 



max 

xe[-2;2]\[a-(5;6+5] 



^s^oo 1/s and therefore large positive 
cr(x), (3.13) 



the Right Hand Side (RHS) of the above equation is well defined provided g[s) G (0; l/mo) 
. However, when g{s) — )• O"*" or when g{s) — )• (l/mo)~, the RHS tends to +oo. Thus, on 
the interval g{s) G (0; \/mo), the RHS must reach a minimum which corresponds to the 
right edge of the density of states. The point g G (0; 1/mo) for which this minimum is 
reached verifies: 



1 



1 



9 



2 + A^„^ 



dx 



[-2;2]\[a-5;b+5] 

and we can compute the right edge of the spectrum s^ 



{l-a{x)gf 



0. 



(3.14) 



1 1 



dx 



[-2;2\\[a-5:b+S] 



,^ from: 

p{x)(j{x) 
1 -a{x)g' 



(3.15) 
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We can now turn to the left edge of the spectrum. Equation (3.12) imphes also the 



asymptotic relation g(s) — )• — c» when s — )• and therefore small positive values of s 



correspond to large negative values of ^(s). The RHS of equation (3.12) is well defined for 
negative values of g{s); it goes to 0~ for very large and negative values of g{s), and goes 
to — oo for g(s) = 0~, so it has a positive maximum somewhere in between. The value of 
this maximum corresponds to the left edge of the density of states and can be computed 
numerically like for the right edge. The point g £ (— oo; 0) for which this maximum is 
reached verifies the same equation as g above, and the left edge Smin is now given by: 

w = i + ^ / dxfi^. (3.16) 

[-2;2]\[a-5;b+5] 

Small fluctuations regime A ^ 5. We first consider the case where A = 6 — a ^ 5, 
corresponding to P ^ Q, in particular the dimension of the perturbed subspace is much 
larger than the dimension of the unperturbed space and so the perturbed space almost 
surely spans the unperturbed subspace. We therefore expect small fluctuations in this 



regime. Equation (3.9) can be solved explicitly in this case. It is in fact possible to perform 



an asymptotic expansion in a{x), which is very small compared to 1 for all x G [—2; 2] \ 



[a — 5]h + S\ and then to solve equation (3.9 ). 



More precisely, in this regime, we have for all x G [—2; 2\\[a — 5]h + S\: 



We plug this approximation in equation (3.9) to obtain 



1 ^ A X p(a) r ^^_ p(x) 



4 ! 



we 



Z{z) N^ J (x - a)2 - A X p{a)Z{z) 

[-2;2]\[a-5-b+S] 

~J^+ / dx-^^^ + Axp{a)Z{z) I dx-^^. 

Z{z) J [x- ay J [x- ap 

[-2;2]\[a-<5;6+<5] [-2;2]\[a-5;b+5] 

Now setting A = J[_2;2]\[a-5-b+S] ^^(i^ and S = A X p{a) /[_2;2]\[a-5;6+5] ^^rfSo^ 
see that Z{z) is solution of the polynomial equation of degree two: 

BZ{zf + {A- z)Z{z) + 1 = 0. (3.17) 

For z = s £ M, this equation has solutions with non-zero imaginary part only if s G 
[A — 2\fB\ A + 2-v/S], which are given by 

-yL + g±V4P-(A-g2) 
^^' 2B 

Using the relation limo^^o QZ{s — iuj) = 7rr(s) for s G M, we find that r{s) in this regime 
is given by the semi-circle law 

r{s) = —— y^iB -{A- s)2, A - 2/b <s <A + 2/B. (3.18) 

2Btt 



This result is consistent with (3.4) since, in this regime, D(a,b;6) = £ A. 

Note that in the particular regime A <^ 6 <^ 1, the quantity B is proportional to A/ 6^ 
and is therefore much smaller than A^ oc 1/5^, meaning that r{s) becoines concentrated 
around s = A, with fluctuations of order y/A/8^. This result is also consistent with the 



EIGENVECTOR DYNAMICS: GENERAL THEORY AND SOME APPLICATIONS 



13 



direct calculation of the root-mean squared fluctuations of s, as obtained in Appendix B, 



see equation (B.2). 




Figure 2. The histogram represents a numerical simulation of the density of 
states of the matrix S (computed with 100 independent samples). The red curve 
is the theoretical corresponding density for r(s) in small fluctuations regime given 
by (3.18). The blue curve represents also the theoretical density r{s) but computed 



numerically by solving directly the system (3.10) and (3.11). For this hgure, we 
chose a = 0,b = 0.01, A = 0.01, S ^ 1. 



Strong fluctuations regime 6 <^ A <^ 1. To simplify notations, we will suppose in the 
following that a and b are such that p{a) ^ p{b). Let us first consider the right edge Smax as 



given by ( 3.15 ) . We need to find an asymptotic expansion in this regime of the g G (0; l/m-o) 



which verifies (3.14). So we start by defining a := gmo G (0; 1) and investigate its behavior 
when 6 <^ A <^ 1. Since mo ^s->-o pio-)/^, equation (3.14) now rewrites as 



p{a)mb 



dx 



p{x)a{x) 



[I 



a 



^M^2 



1. 



(3.19) 



[-2;2]\[a-5;b+5] 

In the limit A ^ 1, it is easy to see that the function a can be written for x < a as 



o[x] 



p{a) 



f 



A 



(3.20) 



where the function / verifies f{u) 



^u^O 



dx 



p{x)a{x) 



(l-«^(^)p{a)^ 



1 and f[u) 



p{a) 



dx 



\/u. Using ( |3.20 ), we can write 



A 



a+2 
A 



du 



[a-x-af{'^)5f 

p{a-uA)f{u) 
{u-ay{u)Y ' 
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where we did the change of variables u = {a—x)/A for the last line. In the limit 5 ^ A <C 1, 
this last integral is dominated by the region where u is small and f{u) ~ 1. We thus have 

ra—5 



(l-aa(x)^)2 A J, {u-aiY 



p{af /■+°° du 
5 1 — a 



Then, using (3.19) and with the same argument now for x > b, we get 

26 



a = 1 



A ■ 



The corresponding g is g = S/p{a){l — 26/ A) and plugging this value of g in (3.15) gives 



P(a) , 1 



dx 



[-2;2]\[a-5;6+<5] 



p{x)a{x) 
1 - a{x)g ■ 



But, it is plain to check that the second term is of order at most ln(J/A)/A ^ 1/6 in the 
limit (5 < A < 1. 
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Figure 3. The histogram represents a numerical simulation of the density of 
states of the matrix S (computed with 20 independent samples) . The red curve 
is the theoretical corresponding density for r(s), it is computed numerically by 



solving the system (3.101 and (3.11). The red dotted vertical lines show the left 



and right edges of the density r{s). The blue dotted curve is the graph of the 
function 8/x'^. For this figure, we chose a = 0, & = 0.1,A = 0.1,(5 = 0.01. 
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Figure 4. The points represent the function s,nax{^) as a function of S. They 
Jly through equations (3.14) and (3.15). 
p{a)/5. In this figure, we chose a = 0, 6 = 0.1, A = 0.1. 



are computed numerically through equations (|3.14p and (|3.15|). The blue dotted 
line is the function 5 



The determination of Smin proceeds similarly, and the calculations are detailed in Ap- 
pendix C. The final result is that 

cp{a) 



where c > is a number of order unity that can be determined if needed. 

To summarize, in this regime, the minimum and maximum eigenvalues Sr, 
of the random matrix E are asymptotically given by: 



and Sr; 



cp{a) 



< s. 



p{a) 



We verified the result for Smax{^) with numerical simulations (see Fig. [4]). 

Together with the exact result on the average value of r{s) in this regime (given by 
D(a,a + A;6)/e'^) and its variance computed in Appendix B, we conjecture that the 
asymptotic behaviour of r{s) in the region Smin "^ s <^ Smax is given by: 



r{s) oc 



(3.21) 



Since the integral of sr{s) is logarithmically divergent (but cut-off at Smin and Smax)-, it is 
easy to see that this form reproduces exactly the logarithmic behavior of D[a, a + A; 5) in 
this regime, see equation (3.5). On the other hand, the integral of s'^r{s) is dominated by 
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its upper bound, leading to a variance of the spectrum given by Smm x Smax, in agreement 



with the exact result obtained in Appendix B, see equation B.l[ Therefore, in this regime, 
the situation is particularly interesting: while most eigenvalues are close to Smin, there is 
a slow power-law tail in r(s) that makes the average of s logarithmically divergent when 
(5 — )• 0. This is why we call this a strong fluctuation regime: the 'overlap' distance D 
between the initial and the target spaces is large because a relatively small number of 
directions are completely lost. 

4. Eigenvector stability for covariance matrices 

4.1. Eigenvectors of Spiked matrices. In this subsection, we will assume that (Cat) 
is a sequence of positive definite matrices. We will denote by A^, . . . , X^ the eigenvalues 
of (Cjv) in decreasing order and we will suppose that 

• there exists a fixed number k < N, q (^ {0; 1) and (Ai > • • • > Afc > (1 + y^)^) 
such that 

[Xi , . . . , Xj.) — s-AT-j.oo (Ai, . . . , Afc). 

• the empirical measure ^n = -^ X]i=fc+i ^>^i converges in the limit of large A^, T 
with N/T = q to the Marchenko-Pastur distribution whose density with respect 
to Lebesgue measure is given by 

1 



Pi^) = 7^ Vil+-^)ix --/-), a<x <b, 

Zirqx 

where 7_ = (1 — ^/q)'^ and 7+ = (1 + y/q)'^. 
For each N, Cat is the true covariance matrix (also called "population covariance matrix"). 
This particular choice for the shape of the matrices Cat is rather natural in view of 
applications. For example, in financial market, the correlation (or covariance) matrix has 
k isolated eigenvalues well separated from the other eigenvalues which form the noisy part 
of the spectrum {Marchenko-Pastur sea or the bulk). 

We now consider the associated empirical covariance matrix E^r defined as: 

1 ^ 



T 

u 



N- 



where the {r\, . . . , r]^), 1 ^ t ^ T are i.i.d. Gaussian vectors of covariance C 

The question we ask in this subsection is: how close are the eigenvectors of E^r to 

those of the matrix Cat ? In the following two paragraphs, we treat the two cases of the 

eigenvectors associated to eigenvalues in the Marchenko-Pastur sea and of those associated 

to the isolated eigenvalues Ai, . . . , A^;. 

This question falls under the scope of section [2] since the matrix Eat can be written as 

a perturbation of the matrix Cat- Indeed we have: 

1 ^ 
^N = Cn + £n, with SN,ij = ^ X] ^^i"^) ~ ^^'*J- ('^■^^ 

t=l 

and the matrix elements of £ are (because of the Central Limit Theorem) of order l/vT 
which is much smaller than 1 as T is large. However, this problem is of different nature 
than the one treated in section [3] because of the non-trivial dependance structure for the 
matrix elements of the perturbation matrix £. It is given by 



£N,ij£N,ke = iCN,ikCN,je + CN,ieCN,jk)/T. (4.2) 
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In the whole current section, • • • denotes an averaging over the r. 



Eigenvectors in the Marchenko-Pastur sea. The results of subsections 3.1, 3.2 and 3.3 
can be extended to this context. We consider the subspace of eigenvectors of Cat cor- 
responding to all the eigenvalues A contained in a certain finite interval [a, b] included 
in the Marchenko-Pastur sea [7_,7+]. We want to compute the distance D between this 
subspace and the subspace spanned by the perturbed eigenvectors of E^r corresponding 
to all eigenvalues of Ejv contained in [a — 6,b + 5], where (5 is a positive parameter. Using 



formula (2.8) as before, we find that in the limit of large N, T (with N/T = q), as soon as 



(5 > 0, the mean overlap distance D is given (using (4.2) for the averaging) by: 



[-2;2]\[a-<5;fe+5] 

where p(A) is now the Marchenko-Pastur distribution of parameter q = N/T. Obviously, 
when T is infinite, D = since E^r = Cm- 

For the singular value density of states r{s), the resolvent of the matrix 5] defined as 
T, = T{I- GGt) now verifies: 

z- ^ + ^ f dA^^^M^ r44) 

[-2;2]\[a-5;b+5] 

where z/ is defined as i>{X) = X f^ dx ,^_li . As before, it is easy to show that the density 
of states of S is compactly supported and to find numerical evaluations of the left and 
right edges. One can also study the limit shape of the density of states in the two regimes 
A <C 5 and 5 ^ A <C 1, with results very similar to the GOE ones above. 

The matrix GG^^ in this case gives a precise information on the relationship between 
the eigenvectors of the population covariance matrix (or true covariance matrix) Cn and 
the eigenvectors of the sample covariance matrix En . Previous works along these lines can 
be found in ^\5\. 

Isolated eigenvectors. In this paragraph, we now consider the case of eigenvectors associ- 
ated to isolated eigenvalues Ai, . . . , A^. We denote by \4>i), \4'2), ■ ■ ■ , \4'k) the corresponding 
eigenvectors of Cat and by \(l)\) , \4>2) , . ■ ■ , \4>'i.) Plthe corresponding eigenvectors of E^v. 

To understand precisely how the \(j)'^) decompose on the basis of the \(j)j) in the limit of 
large N, we want to compute the limit of the average local density of states for each state 
\(p'^) (1 ^ i ^ k), that is the probability measure 



1^ 






(^)^^>^('^'.l'^^-)''^(^-^.-) 



where ^^^ denotes an average over E^v- This expresses the way \(j)'-) is scattered over the 
unperturbed eigenvectors. 

Perturbation theory again allows to compute the quantities {(p'il'pj)'^ for i y^ j: 

[9im) ^^N _ ^Ny T (Af - Af )2 ^^^^ T (Ai - A,-)2 ' 



The vectors \4>i) and \4>'i) depend on A'^ but to simplify notations, we drop the subscript. 
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and for i = j, 








~ Z^ (\N \N\2 ^A'^oo '- 


1 ^ AjAfc 


i'P'im' -- 



Note that the random variables {<j)i\£\(l)j),i ^ j are uncorrelated. Thus, the local density 

of states v^ has k atoms and (for large A^, T with N/T = q) admits a continuous density 
in the Marchenko-Pastur sea. The atom are localized on the \j,j = l,...,k and have 
weights ^ (xlx)2 ^r j ^ i. The continuous density in the Marchenko-Pastur sea [7^,7+] 
is given by: 

j^ .^,^^p P(A)rfA, 7_<A<7+. (4.5) 

(i) 

This asymptotic for the probability measure uj^ has been verified with numerical simula- 
tions. 

4.2. Stability of eigenspaces. We now want to characterize the stability of the subspace 
spanned by the eigenvectors associated to the (largest) isolated eigenvalues. The theory we 
develop here provides a precise estimate of the amount of eigenspace instability induced by 
measurement noise. This sets a benchmark that will allow us to detect any extra dynamics 
of the eigenvectors of the correlation matrix of stock returns in financial markets not 
explained by measurement noise and therefore attributable to a genuine evolution of the 
market (see section [6|) . 



As shown by Eq. (4.1 ) above, the sample covariance matrix Bjis a perturbed version of 
C. Using again the framework of sectional one can calculate the distance (or overlap) be- 
tween the top P eigenvectors of the true correlation matrix C and the top QQ eigenvectors 
of the empirical correlation matrix E. 

Provided T is large enough for the above perturbation theory to be valid, and upon 
averaging over the measurement noise, one gets the following expression for the overlap 
distance D: 

i=lj=Q+l ^ ' 1' 

where the AjS are the eigenvalues of C, in decreasing order. 

Note that one can extend the previous result ( |4.6[ ) to the case where the vectors 
{r\, . . . ,r\j),t ^ are distributed according to a multivariate Student distribution with 

degrees of freedom and covariance matrix C. In this casq [ Eq. (4.6) becomes 



V- 



Note that the Gaussian case corresponds to i^ — t- 00. For v — )• 4^, on the other hand, 
fluctuations become divergent. 

In practice for applications (see sectionp|, one does not know the true correlation matrix 
C and thus it is in fact not possible to compute empirically the overlap distance between 



As A'^ does not have to be necessarily large in this subsection and in the next section, we drop the 

subscript N for the matrices C and E. 

We take Q ^ P as before in section 2 . 
in ^ „ . — ^ . LJ 



10, 



see e.g. Eq. (9.28) p. 154 of [13] that replaces Eq. (4.2 I above 
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the eigenvectors of C and the eigenvectors of the empirical correlation matrix E. However, 
if one is given a time series of empirical correlation matrix (E*)^ ^ o defined for all t as 

1 ^ 

«=1 

where, (rj", . . . , r]^), v ^ are independent Gaussian vectors of covariance matrix C, one 
can similarly define the distance between the eigenspaces of two independent sample co- 
variance matrices E* and E* (determined on two non overlapping time periods, i.e. such 



that |i — s| > T). In this case, the above formula Eq. (4.6) is simply multiplied by a factor 
2. 

For the comparison between the eigenvalues of E* and E*, one can show using per- 



turbation theory (see equation (2.2) and also equation (4.2) for the averaging) that the 



measurement noise is, for T large enough, given by: 

4\2 

(A|-A*)Vs|>T«^- (4.8) 

where the Aj are the eigenvalues of the matrix C measured empirically using the whole 
period of time and where ^^^u_s|>t denotes an empirical average over all s,t such that 
|t — s| > T. As before, if the vectors {r\, . . . ,r'^),v ^ are distributed according to a 
multivariate Student distribution with z/-degrees of freedom and covariance matrix C, one 



finds an extra multiplicative term {u — 2)/{v — 4) in (4.8). 

Another characterization of the stability of eigenspaces was proposed by Zumbach |16j . 
The idea here is to study the stability of the spectral projectors associated to the top k 
eigenvalues. The spectral projector of rank k associated to the top k eigenvalues is defined 
as follows: 

k 
Xk = ^\(t>i){(t)i\ , 

i=l 

where the |0j), i £ {1, . . . ,k} are the eigenvectors of C. As before the true spectral projec- 
tor Xk is measured through an empirical covariance matrix E and the resulting spectral 
projector x'k "^ill be affected by measurement noise. The aim is again to compute proper- 
ties of this spectral projector Xk, so as to be able to separate the measurement noise effect 
from a true temporal evolution of the matrix C. 



Using perturbation theory in Eq. (4.1), we have: 

1=1 i=l \ jf^i -^ 

1=1 J7^« 



+EE«M(i0.)(^.i + i^.)(^.i)+EEE (t'-'Mlt'-'t1 '^^-^^^^' 

i=l j/j i=l jj^i £j^i ^ * J/V « 



where 



1 I ^ {(t>,\£m{^e\£\(j)i) {(^i\£\^i){(^i\£\^,) 



,^^ Y 

Xi - \j \jr' -^i - A£ Xi - Xj 



if^i 
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Using again equation (4.2) 



t>j\£\(t)i) 



we get: 



4 = E i-?E 



l£l 



A^ A j 



if ^ ^ j, 

XjXi/T iij = ij^i^ 
2A?/T otherwise, 



4 = 1 






«)(</'*i + tEE 



1=1 j^i -" 



^'^' m{^,\- 



We see that the vectors (j)i,i £ {1, . . . , N} are also eigenvectors of x^, but with shifted 
eigenvalues. More precisely, we have, for i ^ k 



X'k\^^) 



1 ^ 
T 2-^ 



XiXj 



T ^ (Ai-A,)2 
j=fc+i ^ ^ -I' 



and, for i > k, 



x'k\(Pi 



XiXj 



T 



^,i^^-^j)' 



(4.9) 



(4.10) 



Therefore, in the absence of measurement noise (i.e. for T — )• cx)), x'k ^as k eigenvalues 
exactly equal to unity, and N — k eigenvalues equal to zero, as expected since in this 
case x'k ~ Xk- AH the above results will be compared with empirical data (for the case of 
financial markets) in section 6.1 below. 



5. The case of an isolated top eigenvalue 

5.1. A Langevin equation for the top eigenvalue and eigenvector. A more detailed 
characterization of the dynamics of the top eigenvalue and eigenvector can be given in the 
case where this top eigenvalue is well separated from all the others, as is well known 
to be the case for financial covariance matrices. The financial interpretation of this large 
eigenvalue is the so-called 'market mode': in a first approximation, all stocks move together, 
up or down. In this subsection, we assume that the true covariance matrix C has one large 
eigenvalue Ai of order N well separated from the other ones, which are all equal to A2. We 
suppose that Ai S> A2. 

Let (rf), ^.^., ,l^t^rbe i.i.d. Gaussian vectors of covariance C. Both for technical 

\^/1^2^A'' 

convenience and to follow market practice, we suppose that the covariance matrix is now 
measured through an exponential moving average of the r*. This means that the matrix 
E evolves in time as: 

Eij-,t = (1 - e)Eij-i_i + ery^. (5.1) 

We address the following question: what is the dynamics of the top eigenvalue Ai(i) 
and of the top eigenvector (p\ of the empirical covariance matrix Ej? Of course, the largest 
eigenvalue and eigenvector of the empirical covariance matrix will be, as discussed at length 
above, affected by measurement noise. Can one make predictions about the fiuctuations 
of both the largest eigenvalue and the corresponding eigenvector induced by measurement 
noise? We shall see that such a decomposition is indeed possible in the limit where Ai » A2. 
The calculations in this section and in Appendix D follow closely those made in [T7j which 
were slightly incorrect (see below). 
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We keep the same notations as in the previous section for the eigenvalues of C. The 

li • • • )<f'Ar- 



eigenvalues and eigenvectors of Ej will be respectively denoted as A* , ... , A^ and ^^ ^^ 



Standard perturbation theory, valid for e <C 1, gives: 

A* = {l-e)x\-^ + e{ct>\-'\C\<t>\-') + e{cf>\-'\r^t\^\~'), 

with rjij = TiVj — Cij. Because the returns are Gaussian, we have: 

VijVke = CikCji + CiiCjk- 

In the limit where Ai becomes much larger than all other eigenvalues, the above equation 
simplifies to: 

A* « (1 - e)X\~' +ecos^{et-i)Xi [1 + ^t] , (5.2) 

where cos{9t) = (0* |</>i) and ^t is a random noise term of mean zero and variance equal to 
2. In the limit of large matrices and e — )• 0, the above difference equation can be written 
as a Langevin (or stochastic differential) equation, in the ltd sense: 

dA*i = e \{\i - cos'^{et)\\)dt + V2A1 cos2(6't) dBt] . (5.3) 

where Bt is a standard Brownian motion. We have neglected in the above equation a 
deterministic term equal to esin^(0t)A2, which will turn out to be a factor A2/A1 smaller 



than the terms retained in Eq. (5.3). As we shall show below, the angle 6t turns out to 
be small, so that one can replace cos{6t) by unity in the above equation, which becomes a 
simple Ornstein-Uhlenbeck process. We therefore find for the variogram of Ai: 

' (A^ - A*)') « 28X1 (1 - exp(-e|t - s\)) , 

a result that we mentioned in the above section 14.21 



A similar SDE can be written for the projection of the instantaneous eigenvector \(j)\) on 
the true eigenvector |(/>i). This can again be done using perturbation theory, as is detailed 
in Appendix D. The quantity cos{9t) is found to be close to 1 when e is small, so we set 
xt = l -cos{9t). 

Keeping only the leading term in the three small parameters e, A2/A1 and xt, we finally 
find the following Langevin equation for xt (in the Ito sense): 



dxt = 2eifi- Xt) dt + e\l2xt{4.xt + y) dBt (5.4) 



with, for N — t- 00, e — t- 0, 



u := --— , with q = eN. 
4Ai 



Equation (5.4) defines a very interesting class of random processes, that we call "Poschl- 
Teller" processes, on which we say more in Appendix E. 

In the continuous time limit, we have therefore established two coupled Langevin equa- 
tions (SDEs) for the top eigenvalue A^ and xt. To leading order and for A^ — )• 00, e — )• 0, 
the stationary solution for the "angle" xt can be computed to be: 



/ 4t 
P(x) oc ' 



JV J^ 

2 / 1 \ 2e 



4x + ^ / \ 4x + 4^ 

^ Ai / \ Ai 



which corrects the result obtained in [T7], and is plotted in Fig.^ From the above Langevin 
equation, it is immediate to see that the average value of x is given by x = ^u. It is nicer 
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to rewrite the stationary distribution in terms of x = x/ ji. The interesting regime is when 
q remains of order unity when N ^ oo and e — )■ 0, in which case: 



P{x) « Ze 



Nf{x) 
2 



fix) 



ln(l + qx) 



+ In 1 + 



1 
qx 



where Z is a normahsation. It is easy to see that f{x) has a minimum for x = 1, or x = ;U 
(corresponding to the most probable value), and that /"(I) = 1/(1 + g). This shows that 
the fluctuations of x around x = 1 are of order ^/{l + q)/N and thus very small in the 
large A^ limit. 



Note finally that according to Eq. (5.3), the largest eigenvalue is on average shifted 
upwards compared to the true value Ai, by a factor « (1 + 2/i) = (1 + |i^). This is the 
analogue of a similar well-known result for flat-window averages of empirical covariance 
matrices - see [181 S! ■ 
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Figure 5. A picture of the stationary probability density P{x) of the process xt 
verifying ( |5.4[ ). The parameters are: e = 1/50, N = 200 (corresponding to q = 4) 
and A2/A1 — 0.02. The vertical blue dotted line shows the position of /i w 0.02 for 
this choice of parameters. 



5.2. Variograms. From the Langevin equation one can easily compute the second mo- 
ment Xj with as initial condition xq = 0. Indeed, using Ito's formula and taking expecta- 
tions, we get: 



xt 



Xq + Ae { fi + e 



■2Ai 



x^ds — 4e(l — 2e) / x^ ds . 
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Computing xj with the same technique, we can solve this ordinary differential equation to 
obtain that 



^2^-4.(1-2.)* ^ 



n{ti 



■ e 



A2 

2Ai 



1 



-Aa- 



(2/i + e^)(xo-M) 



2s 

~2et 



1 



-4e(l-2e)t 



-4e(l-2£)t 



l-4e 

In order to characterize the dynamics of the angle fluctuations, we want to compute the 
variogram of xj, defined as f(r) := {xt+r — xt)"^ for r ^ 0, and in the limit t — t- oo. Using 
the previous computations, we obtain, in the scaling limit: 

9/, . „X /x x2 



v{t) 



4N 



1 



-2eT 



We show in Fig. [6] a numerical simulation of the dynamics of the top eigenvector of a fixed 
matrix C such that A2/A1 = 0.033. The resulting variogram compares very well with the 
above prediction. 
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Figure 6. The plain line represents the function u(r) as a function of r for 
e = 0.002, A^ = 200 {q = 0.4) and A2/A1 = 0.033 = 30. The dotted line is a 
numerical simulation of the semivariogram of Xt in the benchmark case where 
there is a constant in time correlation matrix C. 



However, the above calculation is not particularly useful for financial applications, since 
the "true" top eigenvector |(/>i), needed to define the angle 9t, is in general not known. A 
more appropriate quantity to describe the dynamical fluctuations of \(j)i) is, as suggested 
in [TTj, the function r — )• {(pWc/}-^'^) , which we now study analytically. Let us write \cl)\) as 

\^{)=cos{etMi) + \^i), (5.5) 
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where ip\ is a vector in the eigenspace corresponding to the small eigenvalues A2. Therefore: 

m<\>f^) = COs{0t) COs{0t+r) + {V>W±1 ■ 

Now, it is easy to have an explicit expression for 1^2 by considering the empirical covariance 
(or correlation) matrix E* as a perturbation of the true covariance matrix C, as we did 
above. Standard perturbation theory then gives 



!.i|f*|</).)M ,, .^y (0i|f*| 



2^(A,-A,)^/- fe A,-A, 



where 



+00 

It is clear that the last term of the above expression is exactly \(P2), which enables us to 
obtain: 

(A1-A2) ^ 

But, by noting that: 

T-l 

^^ = (1 - 5)^4 + eY.il- ey (r*+^-^rf --^ - C, 

s=0 



and with the fact that ((/)i|<S*|0i)2 = eAiA2/2, we get that: 



and hence, our final result, to lowest order in /j,: 



{m\^l = (1 - ^t)il - xt+r) + (V'ilv'i+O (5.6) 

^ l-2/i(l-e-"^) . (5.7) 

which is similar to the result obtained in [IT], except that the coefficient fj, was a factor N 



too small in that paper. This result will be compared with empirical data in section 6.2 



5.3. Transverse fluctuations of the top eigenvector. In order to go further and 
describe the evolution of the top eigenvector of E (the so-called "market mode" in the 
context of financial markets), we need to study the statistics of the transverse component 
|(/9^). In order to make sense of the pattern created by these transverse fluctuations, we 
propose to introduce the correlation matrix of the components of | if\_ ) in the eigen-basis 
of the true correlation matrix. We therefore define the following A^ — 1 x A^ — 1 matrix: 

1 ^ 

Fij = -YiipimM^j) {ij ^ 2) 

t=i 
The eigenvalues and eigenvectors of this new correlation matrix (not to be confused with 
the empirical correlation matrix E needed to define |v72)0 ^^^^ entirely characterize the 
transverse fluctuations of the "market mode" . 

In the benchmark case where there is a true correlation matrix C stable in time, one 
can check that: 

1 ^ {<Pi\£'\ct^,) {cl>i\£'\^j) 



^'' r^ Ai-A, Ai-A, 
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What is the eigenvalue spectrum of F for this benchmark case? In our case where for all 
i 7^ 1, Aj = A2, the density of states of this type of random matrix has been studied before 
in the literature (see [l9]). Indeed the random variables {(j)i\£^\4)i) are uncorrelated for 
i ^ j^ their mean is and their variance is given by: 

eAi A 



1^2 



(0i|f*|<A.)2 = ^. (5.8) 

However, the random variables {4>i\£^\(j)i) are correlated in time and thus the density of 
states in the limit of large matrices will not be given by the usual Marchenko-Pastur law. 
Rather, {(j)i\£*\(j)i) follows an auto-regressive linear process, for which the authors of |19] . 
give a precise way to compute the density of states in the limit of large matrices by mean 
of its Stieltjes transform. This probability density depends as expected of the parameter 
N/T but also of the parameter e of the auto-regression. In the case where Aj>i = A2, 
one furthermore expects that the eigenvectors of F are isotropically distributed in the 
A^ — 1 dimensional subspace spanned by |(/>2), . . . I^v)- This means that the transverse 
fluctuations \ip±) of the top eigenvector have no particular structure. 

In the more general context where the Aj for i ^ 1 are not all equal to A2 , the eigenvalue 
spectrum of F must be characterized numerically, see below. 

6. Empirical results 

For the following analysis, we have used the daily returns of several pools of stocks 
belonging to 4 major indices: SP500, Nikkei, DAX &: CAC 40. The number of stocks are 
respectively A^ = 500, 204,30, 39 and the period of interest is 2000 — 2010 (11 years of data, 
corresponding to ~ 2750 days) . The main issue, as alluded to above, is that the empirical 

determination of correlatioij | matrices requires some measurement time T. If this time is 

too short, the empirical correlation matrix will appear to evolve with time, but this may 
just be due to the measurement noise that one would like to distinguish from a genuine 
evolution of the underlying structure of correlation. If the measurement time is too long, 
on the other hand, one may miss important correlation shifts and get exposed to unwanted 
sources of risk. 

6.1. Stability of eigenspaces. We first determined the empirical variograms ((A| — 
X\)^)\t-s\=T for i = 1, 2, the result (for i = 1) is shown in Figure PH and is found to be much 
larger than the above theoretical prediction, i.e. 4A?/r, shown as a horizontal plain line. 
The fact that the empirical (red) curve starts from for r = and increases to reach the 
stationary noise level at time r = T is simply due to the overlapping between the sliding 
periods. For those figures, we computed the time series of correlation matrices using a 
sliding window of size T = N (recall N is the number of stocks). Thus, for small markets 
like DAX and CAC40, this value is quite small (respectively 30 and 40) and we find that 
the first eigenvalue of the correlation matrix does not evolve too much during the following 
(non overlapping) period r G [T; 250] days. After this time period, the evolution appears 
and from this point, the difference between the two non overlapping periods increases 
significantly with the time lag. For larger markets such as SPX and Nikkei, the value of 
T is quite large as A^ is respectively equal to 500 and 200. So it is not very surprising 
the temporal evolution shows up immediately. This clearly shows that there is a genuine 
evolution of the eigenvalues of C with time. For the top eigenvalue, this is a well known 



Our results in the previous sections hold for empirical covariance matrices. Hence we centered and 
normalized our empirical time series of returns so as to use them. 
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effect (see |17j and section 5.2 Fig. p^ below): both the volatihty of individual stocks and 
the average correlation between stocks are indeed time dependent, and tend to increase in 
crisis periods |1H llOj. We see that the same is true for smaller eigenvalues too, reflecting 
the instability of intra-sector correlations (data not shown). 
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Figure 7. Plot of ((Af — \\)'^)\t-s\=T as a function of r for the four different 
indexes of our sample. The empirical correlation matrices are computed on a 
sliding window of size T = N . The red line corresponds to the empirical datas from 
our pools of stocks, the plain blue line is the theoretical prediction 4A^/r (valid in 
the limit of large T) and the dotted blue line represents a numerical simulation of 
the benchmark case. Very similar curves hold for the second and third eigenvalues 
as well. 
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But what about the eigenvectors? One could be in a "mixed" situation where the eigen- 
vectors of the true underlying covariance C keep a fixed direction through timq [while its 

eigenvalues are moving around. But if the eigenvalues of the matrix C (that was always 
supposed not to depend of time in the previous sections) themselves are evolving with 
time, the formulas derived in the theoretical section above need to be upgraded. Let us 
assume that the true covariance matrix Cj has time dependent eigenvalues A^ , . . . , A^ 
but with constant eigenvectors that will be denoted |0i), . . . , |0Ar) as above. For times 
s < t with |t — s| ^ T, we define the overlap matrix G**'* as: G^'- = {cj}f\cj}'j) . Under the 
assumption that the eigenvalues are varying sufficiently slowly with time, one now finds 
that: 

D{P, Q;s,t) = -^ (in I det(G^'*^G^'*) 



Up to corrections of order T"^", one can replace in the above formulas the A'^'* by their 
empirical estimates. We finally compute the theoretical distance Dth{P, Q,t) as an average 
over all s,t such that |t — s| = r of the above quantity. 

We now compare our null hypothesis formula, Eq. ( |6.1| ) with (a) an empirical de- 
termination of Demp{P,Q,T) using financial data and (b) a numerical determination of 
Dnum{P, Q, t) using synthetic time series of returns that abide to the hypothesis of a co- 
variance matrix C( with fixed eigenvectors, but time dependent eigenvalues. To achieve 
this, we choose an arbitrary (but fixed) set of orthonormal vectors | ■01 ),..., IV'Af) and 
define Ct as C^ = X]j=i '^ilV'«)(V'i|) where the A* are the empirical eigenvalues obtained 
on the financial return time series. We then use Ct to generate synthetic Gaussian mul- 
tivariate returns {ri('u)}. We show the corresponding results in Fig. ^ with the choice 
P = 5,Q = 10, as a function of r and for T = N days. As above, the study concerns the 
same 4 different pools of stocks corresponding to 4 major indices: SP500, Nikkei, DAX, 



CAC 40. We conclude that (i) the theoretical formula Eq. (6.1) is indeed in very good 
agreement with the numerical results obtained with synthetic data: Dnum ^ Df^; whereas 
(ii) the financial data clearly departs from the null hypothesis of constant eigenvectors, 
since D^mp > Dt^- The same conclusion holds for different values of P,Q. 

We have also computed the value Demp{T = T) for different values of T for every pool 
of stocks, the result is shown in Fig. [9| We compare the empirical function T — )• Demp{T) 
with the theoretical value Dth{T) in the benchmark case where the stock returns are 
distributed as Gaussian vectors of constant covariance matrix C. At first sight, the noise 
contribution appears to be too small to explain the value of Demp{T) at small Ts, at least 
for the pool of the CAC40 and DAX indices. Nevertheless, if we now compare the value 
of DempiT = T) for small value of T with the value of Dth{T = T) in the benchmark 
case where the stock returns are distributed with a multivariate Student distribution with 
i/-degrees of freedom and with a constant covariance matrix C, we see that we can make 
the two curves coincide for small values of T. Therefore, the initial decline as T increases 
indeed follows from a reduction of the measurement noise. However, when T becomes 
very large, the "true" evolution of the eigenvectors starts being visible, and leads to an 



Here we mean that the non-perturbed (or population) eigenvectors do not evolve with time; obviously 
we do not talk about the sample eigenvectors of the empirical covariance matrix E which will be affected 
by measurement noise, evolving around the population eigenvectors. 
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increase of D^mp- This plot suggests that the optimal time scale to measure the empirical 
eigenspaces is around T* = 600 days for the stocks from the Nikkei index, T* = 400 days 
for the ones from CAC40, T* = 450 days for the ones from DAX and T* = 700 days for 
the ones from the SP500 index. 
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Figure 8. Plot of Dth,DnuTn,Demp ior T = N, P ^ 5,Q = 10 for the four 
indices considered here. The blue lines are theoretical benchmark results for fixed 
eigenvector directions (plain line: analytical result, dotted line: numerical simula- 
tions, while the red line is the empirical result). These plots clearly show that the 
subspace spanned by the 5 top eigenvectors evolve with time. 



The above results are fully confirmed, and made more precise, by the spectral projector 



analysis proposed by Zumbach. In Fig. 10 we plot, as in [16], the eigenvalues of the average 



spectral projector x'k ^s a function of its theoretical rank k, for several values of k. We show 
in plain lines the eigenvalues of the empirically determined x'k f°^ ^^^ Nikkei idex, where 
the averaging is made over (overlapping) periods of length T = 600 days, and in dotted 



lines the corresponding theoretical predictions Eqs. (4.9) and (4.10) for the benchmark 
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Figure 9. Plot of D^mpir = T) (red line) and Ah(T = T) as a function of 
T, P = 5, Q = 10 for the four indices considered here. The dotted green line 
represents Dth [t = T) in the benchmark case where the returns are Gaussian with 
constant covariance matrix C and the dotted blue line represents Dth{T = T) in 
the benchmark case where the returns are distributed with a multivariate Student 
distribution with v-degrees of freedom and with a constant covariance matrix C. 
The constant v is chosen equal to 5.5 for the CAC40 and DAX indexes and to 18 
for the Nikkei index. The initial decline as T increases follows from reducing the 
measurement noise. However, when T becomes very large, the "true" evolution of 
the eigenvectors is being felt, and leads to an increase of Demp- This plot suggests 
that the optimal time scale to measure the empirical eigenspaces is around two 
years (T* = 500 days). 



case where the eigenspaces are fixed in time, but are blurred by measurement noise. Here 
again we find clear signals of a true evolution of the eigenspaces. The results for other 
stock indices are very similar. 
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Nikkei 




Figure 10. The dotted lines represent the eigenvahies e^ of the spectral projector 
of rank k (for k — 5, 10, 20 see the legend) as a function of \og{i),i — 1, . . . ,M 
in the benchmark case where the true correlation matrix C is not evolving but 
dressed by measurement noise. The plain lines represent the same function for 
the empirical data from the Nikkei index (204 stocks between 2000 — 2010). Here 
T = 600. In the ideal case (constant correlation matrix, T — >■ cxd), these functions 
should be step functions: e^ ^ ^ = 1 and Ciyk = 0. 



6.2. The dynamics of the top eigenvector. As explained above, one expects in general 
the top eigenvector to wobble around its "true" direction \4>i). The fluctuations around 
\(pi) have two possible origins: one is measurement noise, the other is the presence of a 
systematic rotation of the top eigenvector due to some financial mechanism. 

As a further check that measurement noise is not enough to explain the observed dy- 
namics of \<J3i), we have studied numerically the average overlap of the top eigenvector 
measured a time r apart: {(j)\\(j)-y^'^). This is an interesting quantity because it does not 
require the knowledge of the true direction {(pi). As shown above, this quantity should 
be approximately given by 1 — 2/x(l — e"^"^) if measurement noise is the only source of 
fluctuations. We show in Fig.[Tl]a comparison between this prediction and empirical data 
on the market mode of the Nikkei index. Here again, we flnd that the decorrelation of 
the top eigenvector is much stronger than the benchmark. The deviation from unity is, 
for r = 350, more than three times larger than the benchmark case, with no signs of 
saturation. 

So there is a genuine motion of the top eigenvector in time. This was already pointed 
out in |llj . where we established empirically that the top eigenvector rotates towards 
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Figure 11. The plain line represents the empirical function {(j){\(t)*^'^) as a func- 
tion of T. The period on which the average is performed has 2336 days starting 
01/01/2000. There are N = 204 stocks from the Nikkei index. The exponential 
moving average is made with a parameter e = 1/50. The true empirical correlation 
matrix C is chosen to be the empirical correlation matrix computed using the data 
on the whole period. For this C, we have Ai « 73 and A2 « 0.7. The beginning 
of the period is used to initialize the exponential moving average. The plain blue 
is a numerical simulation in the benchmark case. The dotted line represents the 
function t — )• 1 — 2/i(l — exp(— er)) which corresponds to the benchmark case when 
there is a constant in time correlation matrix. 



the uniform vector |e) = (1,1,..., 1)/\'N when the market goes down, and away from |e) 
when the market goes up. In order to be more comprehensive and understand in details the 
dominant transverse fluctuations of the top eigenvector, we have studied the correlation 



matrix F defined in subsection 5.3 above. We first determined the eigenvalue spectrum 
of F numerically, both for the benchmark case (with only measurement noise) and for 

From this figure, we conclude that, for the Nikkei index 



12 



real empirical data, see Fig. 
during the period 2000 — 2010, there are 3 (maybe 4) eigenvalues of the empirical matrix 
F that reside outside the spectrum of the corresponding benchmark matrix. This suggests 
that these 3 or 4 modes are real and correspond to true fluctuations of the market mode, 



that contribute to the discrepancy displayed in Fig. 11 above. We are now in a position 
to identify the corresponding eigenvectors, i.e. the directions in which the market mode 
most likely to tilt. 

It is natural to think that these directions should themselves correspond to large eigen- 
vectors of the correlation matrix C. Therefore we look for the decomposition of the top 
three eigenvectors of F (that we call \uJi) , \uj2) , {co^)) in terms of \(j)i),i £ {2,3,4,5}. A 
singular value analysis of the 3x4 overlap matrix shows that that one can indeed explain 
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~ 85% of these three eigenvectors in this way, with: 

\ui) « -0.34 \(t)2) + 0.29 103) + 0.30 |04) + 0.84 l^s) 

\U2) « 0.53 \^2) + 0.45 l^s) + 0.47 1(/)4) - 0.54 l^^) (6.2) 

1^3) ~ 0.77 |(/.2) + 0.40 l^s) + 0.48 \^4) ■ 

This means that all the four top eigenvectors of C contribute to the "tilt motion" of the 
market mode. To check that this result is significant, we ran numerical simulations for this 
singular value decomposition in the benchmark case with a constant correlation matrix C 
chosen as before to be the empirical correlation matrix computed using the whole period 
of time (here the decade 2000 — 2010). The 3x4 singular values analysis now give an 
explanatory power of ss 70%, which is clearly less than the 85% obtained above. Still, a 
large part of this explanatory power seems to trivially come from the non random structure 
of C itself. 

In order to revisit the result found in JlT\, we need to understand the link between the 
uniform vector |e) and the eigenvectors \(j)2), ■ ■ ■ , \4>5) of the correlation matrix C. Thus, 
we look at the orthogonal projection |e_L) := (|e) — {e\(pi)\(pi))/M (A/" is chosen such that 
(ej_|ej_)^ = 1 ) of the uniform vector |e) in the space generated by the \4>i),i ^ 2. The 



overlap {e±\(f>i) for all z > 2 are shown in fig. 13 for the Nikkei index during the period 
2000 — 2010. We see that \e±) has indeed very strong overlap with \(p2), I'i^'s); \4'4)j \4'5)j and 
hence, from the above results, also with \uJi), |i^2), l^^s)- Therefore, the fact that the main 
fluctuation modes of {(pi) are along these three uj directions is compatible with the tilt 
motion towards |e). However, other modes, not mentionned in (11) . are detected by the 
present analysis. 

7. Conclusion & Open problems 

Let us try to summarize what we have achieved in this paper. We have developed general 
tools to describe the dynamics of eigenvectors under the influence of small random pertur- 
bations and to study the stability of the subspace spanned by P consecutive eigenvectors of 
a generic symmetric matrix. This problem is relevant in various contexts, including quan- 
tum dissipation and financial risk control, but hopefully the ideas and methods introduced 
here can be used in a much broader context. 

We argue that the problem can be formulated in terms of the singular values of the 
overlap matrix between the initial eigenspace and the target eigenspace, that allows one 
to define an overlap distance, which is small if most of the initial information is conserved. 
We first specialize our results for the case of a Gaussian Orthogonal Ensemble, for which 
the full spectrum of singular values can be explicitly computed in the limit of large matrices 
under the regime where the entries of the perturbation are very small compared to the 
mean level spacing of the non-perturbed matrix. We argue that our setting with rectangular 
Q X P overlap matrices G allows to extend our results to perturbations with entries larger 
than the mean level spacing. We provide some numerical evidences that it is indeed true. 
We find two regimes, depending on the dimension of the target space Q compared to that 
of the initial space P. If Q ^ P, all singular values are close to one another, and their 
distribution is given by Wigner's semi-circle. If on the other hand (Q — P)/P <^ 1, the 
singular values s are distributed according to a very broad law that decays as s~^. These 
results are actually universal, and apply for other matrix ensemble as well - for example 
the case of empirical covariance matrices - provided one is interested in eigenspaces deep 
in the bulk. 
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Figure 1 2 . The red curve represents the cumulative distribution of the density 
of states of the matrix F for the Nikkei index with N — 204 stocks, in the period 
2000 — 2010, with e = 1/50. The blue curve is a numerical simulation for the 
benchmark case with the true correlation matrix C chosen to be the empirical 
correlation matrix using the whole period. For this period and pool of stocks, we 
have Ai ~ 73 and A2 ~ 0.7. 



We have also studied the case of isolated eigenvalues, that are usually very important for 
applications, for example in finance. In most cases, empirical correlation matrices are noisy 
measurements of the true covariance matrix, that can lead to an apparent evolution of the 
top eigenspace, whereas in reality the underlying process is stationary. We have derived 
exact expressions both for the overlap distance and for the average spectral projectors (in- 
troduced by Zumbach [16]) that can be directly compared to empirical results. The special 
case where the top eigenvalue is much larger than all the other ones can be investigated 
in full detail. In particular, the dynamics of the angle made by the top eigenvector and 
its true direction defines an interesting new class of random processes, for which we have 
provided explicit analytical results. 

When compared to empirical correlation matrices of several major stock markets, our 
results allow us to unambiguously conclude that there is a genuine evolution in time of 
the true underlying correlation matrix: measurement noise in itself is unable to explain 
the observed variability (in time) of the top eigenspaces. We have found that the overlap 
distance is minimized when the measurement time is on the order of two to three years. 
Both for shorter and longer averaging times, measurement noise and the genuine evolution 
of the market leads to an instability of the correlation matrix, and to exposure to unwanted 
sources of risk. 

The case of the top eigenvector of the correlation matrix, usually called the market 
mode, is particularly interesting. We have suggested a characterization of the evolution 
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Figure 13. Plot of the overlap \{e±^\4>i) 
shows that the main contribution to |ej^) 
correlation matrix C. 



as a function of i for i ^ 2. This graph 
comes from the top eigenvectors of the 



of its direction through a new correlation matrix, which measures the amplitude of its 
fluctuations transverse to its average direction. We found that the dominant modes are in 
the space spanned by the largest eigenvectors of the correlation matrix itself. 

Now the genuine evolution of the correlation structure of stock returns is well charac- 
terized, one should aim at devising quantitative models for this evolution. As usual, there 
are two ways to do this. One is to postulate an econometric model and try to calibrate it 
on data. In this line of thought, extensions of the GARCH framework have been proposed: 
multivariate GARCH, BEKK model, etc. [20], but they often lack intuition (to say the 
least) and are very hard to calibrate (the a priori number of parameters is of order A^^!). 

The second approach is to think about mechanisms that can lead to changes of the 
correlation structure. For example, market drops may lead to panic sell-offs, that increase 
the top eigenvalue of the correlation matrix and tilt the top eigenvector towards uniformity, 
as reported in |10^ lllj . The impact of rebalancing or deleveraging complex portfolios can 
also lead to substantial changes in the correlation matrix - see the insightful work of 
Cont and Wagalath [21] in this direction. We hope that the tools provided in this paper 
will help building financially motivated, more efficient models of dynamical correlations 
and, correspondingly, second generation risk models where impact and feedback effects are 
accounted for 1221. 
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Appendix A. Proof of the formula (3.6) 
We need to introduce the two level density of states 

1 ^ 



and to note from equation (2.8) that 



2-P Ja J [-2;2]\[a;b] {^ - X'Y 



DiVo;V^) = '-^ I I ^fA^. (A.l) 



From [6| , we know the asymptotic behavior of the two level density of states in the limit 
of large matrices; more precisely, there exists a function g such that, in the limit of large 

p^iX, A') = giNp{X)\X - X'\)piX)p{X')dXdX' (A.2) 

which is defined as g{r) = ^ — {^ — Jq s{t)dt) s'{r) + s(r)^ with s{r) = '^'"^^^ . One can 
check that: 

2 

• in the neighborhood of 0, g{r) ~ ^r, 

• g(r) tends to 1 when r goes to oo, 

• 9'{^) = 0(l/r^) in the neighborhood of oo. 

We can write: 

nr /. A n^ ^^' r f giNp{x)\x-y\) ....,, 
D{a, b;6 = 0) = —— / / ^ p{x)p{y)dxdy. 

^P J a J[~2;2\\[a;b] [^ - VY 

We want to do an asymptotic expansion of the right hand side when N ^ oo. 
First, note that N/P tends to 1/ / p. 
For the integral, we begin by doing an integration by part, we get for x G [—2; a]: 

'' g{Np{x)\x-y\) ^ p{a)g{Np{x){a - x)) p{b)g{Np{x){b - x)) 
[x — y)'^ a — X — X 

+ / ^^ [Piy)9 {Np{x){y - x)) + Np{x)p{y)g' {Np{x){y - x))] . (A.4) 



We need to integrate equation (A. 3) between —2 and a and between b and 2 and to 
compute the asymptotic of every integrals of the right hand side. We will decompose each 
integral into two terms so as to take advantages of the asymptotic property of g around 
and oo. 

Set rj = N~^~^" with a > 0. First we consider the integral: 

, .g(Np(x)(a — x)) , f'^" dx , x. , , x. . 



/■^" dx 
pi"-) / —9ipia)x) = p{a) 
Jo a; 



ln(p(a)A^")- / ln{x)g'{x)dx 
Jo 
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Using the fact that g{r) tends to 1 when r goes to oo, we easily get that, in the Hmit 

N -^ oo: 

-2 a -a; 

~ / ■^^dx = -p{a)ln{N-^+")+ f p'{x)ln{a-x)dx. 



-2 a- X 
Moreover, we easily find that, when N tends to oo: 

-2 b-x 



" p{x) 



b — X 



dx . 



and 



, ^giNpix)ib-x))^ ^ r p{x) ^ 
p{x) ^, ^\ -dx^ / T^^dx, 



b — X 



b — X 



which goes to as 77 goes to 0. 

The next term is easy to control using the fact that g{r) goes to 1 when r goes to 00; 

as A^ —7- 00: 

dxp{x) g{Np{x){y-x)) -' -J~.-~^~.^ 

-2 Ja y-x 



dxpix) 

2 Ja y-x 



Using the fact that g'{r) is of order 1/r^ for large r, it is easy to check that 

/"""'' /"^ dy 
N dxp{x) p{y)g' {N p{x){y - x)) 

J-2 Ja y ~ X 



[p{y)9 {Np{x){y - x)) + Np{x)p{y)g' {N p{x){y - x))] 



is of order N ". 

The remaining term 

dxp{x) j 

I a—rj J a i/ X 

is of order N^^^"^ . 

One has to go through the same steps to compute the asymptotic of the integrals 
between b and 2. 

Finally, we get: 



D{a,b-5 = Q)^\uN e^ P(«)' + pW% ^( fc) 

2/>(A)dA 



(A.5) 



where 



A{a,b) 



{p{af + p{bf) [I - I \n{x)g'{x)dx 



^laP 

+ p{afHp{a))+p{bfln{p{b)) 



+ p{a) I p{x)ln{a — x)dx — p{b) / p{x)ln{x — b)dx 
-2 Jb 

/, X /" p{x)dx /"^ p{x)dx 



r [^ p' 

+ / dxp{x) I — 

J-2 Ja y 



p'{y)dy 

y -X 



dxpix) I 

J a x-y 
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Appendix B. Derivation of the standard deviation a{r) of r{s). 
We have 

a\r) = {risf) - {r{s)f 



^(tr(S^))-((^tr(S))V 



P 

But the two quantities are computable easily in the limit of large matrices using the 
convergence of the density of states for Hq; We obtain 



1 
P 



e 



(tr(^'))-^ d^ 



dX' 



dX 



„ p{X)piX')piX") 
Nt J ^" J ^" J ^" (A-A02(A-A")2 

[a;b] [-2;2]\[a-5;b+S] [-2;2]\[a-S-b+S] 

e' r .. f .., f .-„ piX)p{X')piX") 



+ ^y'V' J ''\x-x"nx'-x") 

[a;b] [a;b] [-2;2]\[a-<5;6+<5] 



and 



,p(A)p(A')p(A")p(A'") 



(I itr(E))') . -£l / .A / dA' / dA" / dV j^ _ ^„j,^^, _ ^„,j, 

" [a;fe] [a;b] [-2;2]\[a-5;fe+5] [-2;2]\[a-5;fe+5] 



Those two expressions give in the regime 5 ^ A ^ 1 

p{a) 



o-(r) 



/JA 



and in the regime A <ti 6 -^ 1, 



a{r) f^ p{a) 



2A 



3^3- 



(B.l) 



(B.2) 



Appendix C. Determination of Sm,„ in the strong fluctuation limit 



It is given by (3.16) and we have to compute the g G (— oo;0) that verifies (3.14). For 
simplicity, we set g = —Ij and we aim to compute ff ^ such that 



m J 



dx 



p{x)a{xy 



1 



Nl J -{l + a{x)gr 92 ' 

[-2;2]\[a-5;b+5] 

As g is non- negative, the integral on the left hand side converges when S goes to and 
hence g verifies in fact 



dx 



p{x)a{x)^ 



1 



1 

Nb J -(l + a(x)^)2-52' 

[-2;2]\[a;6] 



(C.l) 



We now need to estimate the integral in the limit A ^ 1. As before, we can write using 



(3.20) 



'^ p{x)a{xf _pial f"^ , pia-uA)fiu) 

_2 (1 + a(x)5)2 - A L (..^pMr 



(n+™/(^)) 
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In the limit A <C 1, this integral is dominated by the region where u is small and f{u) ~ 1 
and hence we have the following estimate 



ax — — — ~ 



du 



2 {l + a{x)gY A Jo {u + ^ 

p{af 
9 



)9\2 



Then we deduce from (C.l ) and with the same argument for the integral between b and 2 
that 

._ A 



Now we have to plug this g into equation (3.16) to obtain 

p{x)a{x) 



Ma) I 1 

A ^7V„^ 



dx 



l + a(x] 



2p{a) 



[-2;2]\[a;b] 

To evaluate the integral, we need to cut it into two parts. The first part is handled by 
'■" p{x)a{x) 



2 l + (y{x) 



2p(a) 



"i' p{a - uA)f{u) 



Finally, we can deduce that 



2pia) 

A 



p(a) I du „/ X , 

'^ 'Jo u + f{u)/2 

Jo u + f{u)/2 

f{u) 



+ 00 



du „/ s , 

u + f{u)/2 



1 . 



Appendix D. Derivation of (5.4) 
Using perturbation theory, one gets, in braket notation: 



\^{) 



t-l|y,i^t*Ut-l\2 



'-tE "w:\'^./ i^">^-e 



2 



i^l 



,t-i|^V*|0*-i) 



it-i\ 



e 



'-n^^^:'" 



-i'rV*uri\2 



^1 



iJLl 



Since cos(0t) = {4>\\(j)i), we can write 



(^*i = cos(0i)|(Ai)+sin(0t)|V9i) 
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where \^±) is a vector lying in the subspace spanned by the vectors \4>2)-, ■ ■ ■ ■, I'/'Af)- We 
want to describe the dynamic of cos{Ot)] we deduce from the previous equation that 



e 



2 



cos(^*) = ^1 - ^^3^ {{4>\-\r'r'y\c^\-') - {c^\-'Vr'*\ct>\-'?)j cos(^,-i) (D.l) 
+ T^ ((<Ai|rV*>*f 1) - (<^*f VV*>i-^)cos(ei_i)) . (D.2) 

Since we have: 



(0i|C|(/)*f') = Aicos(0i_i), 



{ct^r\r^tWi-'? = 2cos\dt^i)\i + s\n\dt-i)\i\2 



i"'l^t|0*r')2 = 2 {\^cos\et^i) + \2siYiHet-i)y 



b\-\h^*\<l,\-y = \lcos\et-i) + 2coB\et-i)\l + Bm\et-i)\i\2, 



{^\-Wri*f\(t>\-') = cos^(0t_i)(3Af + {N- l)AiA2) + sin^(^i_i)((iV + 1)A^ + A1A2) , 



equation (D.l) can be rewritten, in the asymptotic regime where e ^ 1,A^ ^ 1 with 
q = eN fixed and A2 <C Ai, keeping up to terms of order 2 for "drift" terms and of order 
1 for noise terms in e and A2/AJ [ 

d{cos{et)) = —^2 [(^? + ^^1^2) cos\et) - \l cos\et)\ cos{et)dt 



+ e cos{9t) sin^ {et)dt + atdBt 



where 



£2 



a2 = -2 [2A? cos^{9t) sm\0t) + A1A2 cos2(26't)] sin2(^t). (D.3) 

When 9t "^ I, this leads to Eq. |5.4| given in the main text for xt = 1 — cos{dt). 



Appendix E. Transition probability of xt 
In this appendix, we show that the function P{x, t) giving the probability that the 



"particle" xt verifying (5.4) is in x at time t can be computed explicitly. More generally, 
we will show that one can compute explicitly this transition density P(x, t) for a process 
Xt with initial condition in t = given by xq ^ verifying the Langevin equation 

dxt = e{ii-xt)dt + ay/j^^{x^+h)dBt (E.l) 

where 0, /i, a and h are positive constants and Bt a standard Brownian motion. One can 
proceed to the change of variables 

yt = cosh~^ (-xt + l\ <^ ^* ^ 2 (^"^'^(y*) ~ -*-) ' 
and find that the process yt verifies 

V h smh(yf) 2 smh(yt) / 



13 



Note that sa?{dt) ~ 2/i is of order A2/A1 and that 1 — cos(^t) ~ ^ is also of order A2/A1 
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We will denote by F{y) the drift coefficient of the previous stochastic differential equation 



(E.2) and denote by U its potential, which verifies U' = —F. The transition density P{y, t) 



verifies the Fokker-Planck equation 

dP _ d{FP) cj2 d^P 
dt dy 2 dy'^ 

By setting P{y,t) := e~^^y>'" il){y,i), this equation becomes a Schrodinger equation: 
with the so-called Poschl- Teller potential V{y): 

1 f a (3 cosh(y) 



2 ^(y)V', 



2 \sinh (y) sinh (y) 



+ 7 



with: 



"H.4)(^^)4(-!^y 



7 = -,{0 + 






/3 = 2«(l + =^)(l + ^ 



Since the evolution of ip{y, t) is governed by a self adjoint operator 

we can use its eigenfunctions to construct an orthonormal basis (V'n) with corresponding 
eigenvalues (— A„). The general solution tl'{y, t) can thus be expanded in the following form 

n 

The general solution for P is thus given by 



The initial conditions for yt determines the sequence (c„). In particular, if at time t = 0, 
the probability P{y, 0) = 6{y — yo) with yo '■= argcosh(|xo + 1), then it is straightforward 
to see that 

The spectrum of Ti consists of a discrete and a continuous branch. The discrete energy 
levels (eigenvalues) are computed in, e.g. [25] and are given for all n G N, n ^ g/2 with 
5 = 1 + 20(72, by 

2 

^n = -Yn{g-n) . (E.3) 

The corresponding eigenvectors are also computed in [25] and are expressed in terms of 
Jacobi polynomials. To the best of our knowledge, the continuous branch of the spectrum 



EIGENVECTOR DYNAMICS: GENERAL THEORY AND SOME APPLICATIONS 41 



has not been fully characterized in the literature. We should also mention that in the limit 
6^0 the corresponding process has been studied in details (see [26j and the appendix of 
[23]). The problem can now be mapped into the Morse potential, which has exactly the 
same discrete spectrum as above (as expected since b does not appear) , with eigenfunctions 
that can be expressed in terms of Laguerre polynomials. However, we have not been able 
to directly match the eigenfunctions in the two cases, and understand the 5 — )• limit in 
details. The limit 6 — )• oo with a'^b fixed, on the other hand, boils down to the standard 
Bessel process with mean-reversion. 
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